Este é um R Notebook. Quando você executa o código R no notebook, os resultados aparecem abaixo do código.
O código R precisa estar em “pedaços” (chunks) em um Notebook R. Logo abaixo, está um exemplo de um pedaço de código R. O código produz uma parábola.
Tente executar este cógdigo, colocando o cursor dentro do código e clicando no botão Run, ou colocando o cursor dentro do código e pressionando Ctrl+Shift+Enter (Windows/Linux).
x <- seq(-1, 1, by = 0.01)
y <- x^2
plot(x, y, type = "l")
Para ocultar a saída, clique nas setas para expandir/recolher a saída. Para limpar os resultados (ou um erro), clique no botão “x”.
Você também pode pressionar Ctrl+Enter (Win/Linux) para executar uma linha de código por vez (em vez de todo o bloco).
Adicione um novo pedaço de código R clicando no botão Insert new code chunk na barra de ferramentas ou pressionando Ctrl+Alt+I (Win/Linux).
Insira um novo pedaço de código R abaixo, digite e execute o código: Sys.time()
Em vez de usar teoria e fórmulas, vamos explorar e revisar a modelaos de regressão linear usando dados simulados.
O código abaixo faz o seguinte:
atribui a x os valores entre 1 e 25.
Em seguida, geramos y como função de x usando a fórmula 10 + 5*x
criamos a data frame d para armazenar os valores de x e y
fazemos um gráfico de dispersão simples entre y e x.
x <- 1:25
y <- 10 + 5*x # formula for a line
d <- data.frame(x, y)
plot(y ~ x, data = d)
Agora vamos adicionar algum “ruído” aos nossos dados adicionando valores (pseudo) aleatórios de uma distribuição Normal com média = 0 e desvio padrão = 10.
A função rnorm() nos permite extrair valores aleatórios
de uma distribuição Normal.
O comando set.seed(1) garante que todos geremos os
mesmos dados “aleatórios”:
set.seed(1)
erro <- rnorm(n = 25, mean = 0, sd = 10)
# Adiciona erro aleatorio a 10 + 5*x e refaz o grafico de dispersao
d$y <- 10 + 5*x + erro
plot(y ~ x, data = d)
Agora y parece estar associado a x, mas não completamente determinado por x.
y é a combinação de uma parte fixa e uma parte aleatória:
10 + 5*xrnorm(n = 25, mean = 0, sd = 10)E se recebêssemos esses dados e nos dissessem para determinar o processo que os gerou? Em outras palavras, trabalhe de trás para frente e preencha os espaços em branco:
Essa é a abordagem tradicional em modelagem/regressão linear. Você tem alguma variável resposta numérica e deseja encontrar o modelo que gerou os dados.
A modelagem de regressão linear múltipla tradicional assume as seguintes hipóteses como verdadeiras (entre outras):
a fórmula é uma soma ponderada de preditores (por exemplo, y = 10 + 5*x).
o erro é uma realização aleatória de uma distribuição Normal com média = 0.
o desvio padrão da distribuição Normal é constante (por exemplo, 10)
Os modelos de regressão linear tentam estimar os “pesos” da primeira hipótese e o desvio padrão da terceira hipótese.
Vamos demonstrar: Abaixo tentamos recuperar o processo gerador para
nossos dados. Para isso usamos a função lm(). Temos que
especificar a fórmula para a primeira hipótese. A segunda e a terceira
hipóteses estão incorporadas em lm().
A fórmula “y ~ x” significa que pensamos que o modelo é “y = intercepto + inclinação*x” ou (“b0 + b1x”).
A menos que especifiquemos o contrário, isso pressupõe que queremos
estimar o intercepto. Isso diz à função lm() para pegar os
dados e encontrar os melhores intercepto e coeficiente de inclinação.
Observe que este é o modelo correto!
Quando você usa lm() você geralmente deseja salvar os
resultados em um objeto. Abaixo salvamos no objeto “mod_sim”. Em
seguida, visualizamos os resultados do modelo usando
summary():
mod_sim <- lm(y ~ x, data = d)
summary(mod_sim)
Call:
lm(formula = y ~ x, data = d)
Residuals:
Min 1Q Median 3Q Max
-23.876 -4.613 2.254 5.909 14.648
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.135 3.999 2.784 0.0105 *
x 5.042 0.269 18.743 1.98e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9.7 on 23 degrees of freedom
Multiple R-squared: 0.9386, Adjusted R-squared: 0.9359
F-statistic: 351.3 on 1 and 23 DF, p-value: 1.981e-15
O modelo retorna as seguintes estimativas:
y = 11.135 + 5.042 * x
erro = rnorm(n = 25, mean = 0, sd = 9.7)
Eles estão bem próximos dos valores “verdadeiros” de 10, 5 e 10 que usamos para simular os dados.
Na vida real, NÃO CONHECEMOS a fórmula da parte 1. O verdadeiro processo gerador dos dados será muito mais complicado. A fórmula que propomos será apenas uma aproximação e pode não ser boa.
Na vida real, NÃO SABEMOS se a suposição de normalidade ou suposição de variância constante do ruído são plausíveis.
Como podemos avaliar nosso modelo?
Poderíamos usar nossas estimativas dos parâmetros do modelo para gerar dados e ver se eles se parecem com nossos dados originais.
Execute o código abaixo de uma vez e mais de uma vez. Os pontos pretos não mudam, mas os vermelhos sim. Isso parece muito bom! Nossos dados gerados pelo modelo parecem semelhantes aos nossos dados observados.
# o modelo simulado original é: d$y <- 10 + 5*x + erro
d$y2 <- 11.135 + 5.042*d$x + rnorm(25, 0, 9.7)
plot(y ~ x, data = d) # dados originais
points(d$x, d$y2, col = "red") # dados simulados
Também podemos comparar curvas de densidade suaves dos dados originais e com os gerados pelo modelo. As curvas de densidade suaves são basicamente versões suaves de histogramas.
Se tivermos um bom modelo, os dados gerados pelo nosso modelo deverão ter uma distribuição semelhante aos dados originais. Execute o código abaixo de uma vez e mais de uma vez.
hist(d$y, freq = FALSE) # freq = FALSE -> area das barras somam 1
lines(density(d$y)) # dados originais
d$y2 <- 11.135 + 5.042*d$x + rnorm(25, 0, 9.7)
lines(density(d$y2), col = "red") # dados simulados
O resultado parece bom. A distribuição dos dados gerados pelo nosso modelo é muito semelhante aos dados observados. Você deve fazer isso mais de uma vez, digamos 50 vezes, para garantir que o modelo gere consistentemente dados semelhantes aos observados. Mostraremos uma maneira mais eficiente de fazer essa simulação mais adiante.
Como achamos que nosso modelo é “bom”, podemos usá-lo para fazer uma previsão. Por exemplo, quando x = 10 qual é o valor esperado de y? Dito de outra forma, qual é a média de y condicional a x = 10?
Podemos obter essa previsão com a função predict(). O
argumento interval = "confidence" signififca que desejamos
uma estimativa por intervalo com 95% de confiança (IC) para esta média
condicional.
predict(mod_sim, newdata = data.frame(x = 10), interval = "confidence")
fit lwr upr
1 61.55932 57.2126 65.90603
A média esperada (condicional) de y quando x = 10 é de cerca de 61,6 com um IC de 95% de (57,2, 65,9). O IC nos dá uma noção de quão incerta é essa média esperada. Na verdade, seria melhor relatar isso como “a média esperada de y quando x = 10 deve estar entre 57 a 66”.
Poderíamos também tentar resumir a relação entre y e x examinando os
coeficientes (ou pesos) obtidos no sumário dos resultados. Podemos
extrair os coeficientes do sumário usando a função
coef():
coef(summary(mod_sim))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.134859 3.9994866 2.784072 1.054874e-02
x 5.042446 0.2690346 18.742742 1.981382e-15
O coeficiente x diz que y aumenta cerca de 5 unidades para cada aumento de uma unidade em x, variando em média 0,27 acima ou abaixo de 5. O erro padrão nos dá uma indicação da incerteza nesta estimativa. Falaremos mais sobre os valores t e valores-p adiante.
**RESUMO:**
Essencialmente, a modelagem de regressão linear básica consistem em:
propor e ajustar um modelo linear
determinar se o modelo é bom e se as premissas são atendidas em sua maioria.
usar o modelo para explicar a relação entre y e x e/ou fazer previsões.
Vamos ver o que acontece quando ajustamos um modelo “ruim”. Abaixo,
adicionamos uma nova coluna à data frame d chamada
z, que é uma amostra aleatória de números no intervalo de
-100 a 100. runif() amostra números de uma distribuição
uniforme entre min e max.
set.seed(4)
d$z <- runif(25, min = -100, max = 100)
LEMBRETE: É possívei adicionar um novo trecho de código R clicando no botão Insert new code chunk na barra de ferramentas ou pressionando Ctrl+Alt+I (Win/Linux).
Modele y como uma função de z usando
lm(y ~ z, data = d) e salve os resultados em um objeto
chamado mod_sim2, veja os resultados usando a função
summary(). Qual modelo foi estimado? Qual é o desvio padrão
estimado do erro aleatório normalmente distribuído?
Use o modelo para simular um histogramas ou densidade e compare
com o histograma e densidade original de d$y. Execute o
código várias vezes para ver como a curva de densidade gerada pelo
modelo varia.
Agora, vamos aplicar um modelo de regressão linear a dados reais.
Vamos importar os dados que usaremos . Os dados com os quais trabalharemos são dados imobiliários do condado de Albemarle no estado da Virgínia/EUA, baixados do Office of Geographic Data Services. Usaremos uma amostra aleatória dos dados.
library(readr)
url <- 'https://raw.githubusercontent.com/clayford/dataviz_with_ggplot2/master/alb_homes.csv'
homes <- readr::read_csv(file = url)
dplyr::glimpse(homes)
Rows: 3,025
Columns: 15
$ yearbuilt <dbl> 1754, 1968, 1754, 1934, 1963, 1754, 1932, 1960, 1950,…
$ finsqft <dbl> 1254, 1192, 881, 480, 720, 990, 792, 1232, 576, 863, …
$ cooling <chr> "No Central Air", "No Central Air", "No Central Air",…
$ bedroom <dbl> 1, 3, 2, 0, 2, 3, 2, 3, 2, 0, 2, 1, 2, 3, 3, 3, 4, 4,…
$ fullbath <dbl> 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 2, 1, 2,…
$ halfbath <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,…
$ lotsize <dbl> 4.9330, 1.0870, 195.9300, 10.0000, 1.0000, 1.0700, 4.…
$ totalvalue <dbl> 124300, 109200, 141600, 69200, 139700, 67400, 237500,…
$ esdistrict <chr> "Brownsville", "Scottsville", "Stony Point", "Crozet"…
$ msdistrict <chr> "Henley", "Walton", "Sutherland", "Henley", "Henley",…
$ hsdistrict <chr> "Western Albemarle", "Monticello", "Albemarle", "West…
$ censustract <dbl> 111.00, 113.01, 104.01, 101.00, 102.02, 112.01, 101.0…
$ age <dbl> 265, 51, 265, 85, 56, 265, 87, 59, 69, 265, 71, 265, …
$ condition <chr> "Substandard", "Substandard", "Substandard", "Substan…
$ fp <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1,…
Vejamos as primeiras 6 linhas:
head(homes)
Dicionário dos Dados:
Digamos que queremos modelar o valor total médio de uma casa
(representado pela variável totalvalue) em função de várias
características, como tamanho do lote, metros quadrados construídos,
presença de ar central, etc.
Vamos ver como a distribuição da variável que representa o valor total usando um histograma. Observe que a distribuição é bastante assimétrica:
hist(homes$totalvalue)
Podemos também criar uma curva de densidade, que é uma versão suavizada de um histograma:
plot(density(homes$totalvalue))
Para modelar o valor médio total das residências em função de diversas características, precisamos propor um modelo linear. Ao contrário do exemplo anterior, estes não são dados simulados para os quais conhecemos o processo gerador dos dados.
Como propor um modelo? É fundamental ter algum conhecimento sobre o fenômeno
Vamos ajustar um modelo linear usando pés quadrados
(finsqft), o número de quartos (bedroom) e
tamanho do lote (lotsize). O sinal de mais (+) significa
“incluir” no modelo.
m1 <- lm(totalvalue ~ finsqft + bedroom + lotsize, data = homes)
summary(m1)
Call:
lm(formula = totalvalue ~ finsqft + bedroom + lotsize, data = homes)
Residuals:
Min 1Q Median 3Q Max
-1164152 -82296 -7690 57164 5879188
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -133328.25 16273.88 -8.193 3.73e-16 ***
finsqft 284.46 5.37 52.967 < 2e-16 ***
bedroom -13218.41 5750.70 -2.299 0.0216 *
lotsize 4268.77 219.32 19.464 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 227200 on 3021 degrees of freedom
Multiple R-squared: 0.6164, Adjusted R-squared: 0.616
F-statistic: 1618 on 3 and 3021 DF, p-value: < 2.2e-16
A função coef() extrai os coeficientes, ou parâmetros
estimados:
coef(m1)
(Intercept) finsqft bedroom lotsize
-133328.2482 284.4613 -13218.4091 4268.7655
Com os parâmetros estimados, podemos escreve o modelo estimado:
totalprice = -133328.2482 + 284.4613*finsqft + -13218.4091*bedroom +
4268.7655*lotsize`
Algumas interpretações básicas:
cada pé quadrado acabado adicional adiciona cerca de US$ 284 ao preço, mantidas as demais variáveis constantes.
cada quarto adicional reduz o preço em US$ 13.218, mantidas as demais variáveis constantes.
cada acre adicional de tamanho de lote adiciona cerca de US$
4.268 ao
preço, mantidas as demais variáveis constantes.
Cada uma dessas interpretações assume que todas as outras variáveis são mantidas constantes! Portanto, estima-se que adicionar um quarto a uma casa, sem aumentar o tamanho do lote ou os metros quadrados acabados da casa, reduza o valor da casa. Isso faz sentido?
Este é um modelo “bom”? Vamos simular os dados do modelo e compará-los com os dados observados. Um “bom” modelo deve gerar dados semelhantes aos dados originais.
Poderíamos fazer isso manualmente:
sim_y <- -133328.2482 + 284.4613*homes$finsqft +
-13218.4091*homes$bedroom + 4268.7655*homes$lotsize +
rnorm(3025, sd = 227200)
Uma maneira mais fácil e rápida é usar a função
simulate() que permite simular múltiplas amostras. Aqui
geramos 50 amostras. Cada amostra terá o mesmo número de observações que
nossa amostra original (n = 3.025). Cada valor de amostra é gerado
usando nossos valores observados para finsqft,
bedroom e lotsize. O resultado é uma data
frame com 50 colunas.
sim1 <- simulate(m1, nsim = 50)
Agora vamos representar graficamente os dados simulados e os dados
observados usando gráficos de densidade. Usamos um loop for
para adicionar estimativas de densidade suaves das 50 simulações.
O comando sim1[[i]] extrai a coluna i como um
vetor. (execute o código de uma vez.)
plot(density(homes$totalvalue))
for(i in 1:50)lines(density(sim1[[i]]), col = "grey80")
(Veja o final do notebook para saber como criar esse gráfico usando o pacote ggplot2 e como transformar esse código em uma função.)
Este não parece ser um bom modelo. Na verdade alguns dos nossos valores simulados são negativos!
Antes de revisarmos o modelo, lembre-se das principais hipóteses:
totalvalue pode ser modelado por uma soma ponderada:
valor total = interceptação + pés quadrados + quartos + tamanho do
loteA linguagem R fornece alguns gráficos de diagnóstico básicos para
avaliar as hipóteses 2 e 3 sobre os resíduos. Basta aplicar a função
plot no objeto que armazenou os resultados da estimação do
modelo.
plot(m1)
Como interpretar os gráficos:
Resíduos vs Valores Ajustados: deve possuir uma linha horizontal com dispersão uniforme e simétrica dos pontos; se não, haverá evidência de que a variância ou desvio-padrão não é constante.
QQ normal: os pontos devem ficar próximos à linha diagonal; caso contrário, haverá evidência de que o resíduo não segue, aproximadamente, uma distribuição n=Normal.
Scale-Location: deve ter uma linha horizontal com dispersão uniforme de pontos; (semelhante ao nº 1, mas mais fácil de detectar tendência na dispersão).
Residuais vs Alavancagem: pontos fora das curvas de nível são observações influentes. Alavancagem é a distância do centro de todos os preditores. Uma observação com alta alavancagem tem influência substancial no valor ajustado.
Por padrão, os 3 pontos “mais extremos” são rotulados pelo número da linha. 2658 aparece em todos as quatro gráficos. É uma casa muito grande e cara.
homes[2658,]
Esses gráficos revelam que nossas suposições sobre a normalidade e variância constante dos resíduos são altamente suspeitas. Nosso modelo é simplesmente ruim.
O que podemos fazer?
A variação não constante pode ser evidência de um modelo errado ou de uma variável resposta muito assimétrica (ou um pouco de ambos). Lembre-se de que a variável resposta é bastante assimétrica:
hist(homes$totalvalue)
Ao lidar com uma variável resposta estritamente positiva e muito assimétrica Ao lidar com uma resposta estritamente positiva e muito distorcida (como variáveis que representam preços), é comum transformar a variável resposta para uma escala diferente.
Uma transformação comum é uma transformação logarítmica. Quando
aplicamos o logaritmo natural a variável totalvalue, a
distribuição parece um pouco mais simétrica, embora seja importante
observar que isso não é uma suposição da modelagem linear!
hist(log(homes$totalvalue))
Vamos tentar modelar a variável totalvalue
log-transformada.
m2 <- lm(log(totalvalue) ~ finsqft + bedroom + lotsize, data = homes)
Os gráficos de diagnóstico dos resóduos parecem melhores.
plot(m2)
Mas este é um “bom modelo”? Nosso modelo proposto de somas ponderadas é bom? Novamente, vamos simular dados e comparar com os dados observados.
sim2 <- simulate(m2, nsim = 50)
plot(density(log(homes$totalvalue)))
for(i in 1:50)lines(density(sim2[[i]]), lty = 2, col = "grey80")
Os resultados não são tão ruins!
Digamos que estejamos satisfeitos com este modelo. Como interpretamos as estimativas dos parâmetros? Como a variável resposta foi transformada usando a função logaritmo natural, interpretamos as estimativas dos parâmetros como diferenças proporcionais aproximadas.
Abaixo vemos os coeficientes arredondados para 4 casas decimais.
round(coef(m2), 4)
(Intercept) finsqft bedroom lotsize
11.6189 0.0005 0.0429 0.0047
Estas são proporções. Para obter porcentagens, multiplique por 100.
round(coef(m2), 4) * 100
(Intercept) finsqft bedroom lotsize
1161.89 0.05 4.29 0.47
Algumas interpretações básicas:
cada metro quadrado construído adicional aumenta o preço, em média, em 0,05%. Ou multiplique por 100 para dizer que cada 100 pés quadrados acabados adicionais aumenta o preço em 5%.
cada quarto adicional aumenta o preço, em média, em cerca de 4,3%.
cada acre adicional de tamanho de lote aumenta o preço, em média, em cerca de 0,47%.
Lembre-se, a interpretação assume que todas as outras variáveis são mantidas constantes!
Uma estimativa um pouco mais precisa pode ser obtida
exponencializando os coeficientes e depois interpretando os
efeitos como multiplicativos em vez de aditivos. Abaixo
exponenciamos usando a função exp e depois arredondamos
para 4 casas decimais.
round(exp(coef(m2)), 4)
(Intercept) finsqft bedroom lotsize
111177.8502 1.0005 1.0439 1.0047
Por exemplo, cada quarto adicional (assumindo que todas as demais variáveis preditoras se mantém constante) aumenta o preço total esperado em cerca de 4,4%. Multiplicar por 1,0439 equivale a somar 4,39%.
Vamos revisar o resultado do resumo:
summary(m2)
Call:
lm(formula = log(totalvalue) ~ finsqft + bedroom + lotsize, data = homes)
Residuals:
Min 1Q Median 3Q Max
-2.94133 -0.14900 0.02095 0.16882 2.60752
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.162e+01 2.389e-02 486.272 < 2e-16 ***
finsqft 4.764e-04 7.885e-06 60.413 < 2e-16 ***
bedroom 4.292e-02 8.443e-03 5.083 3.94e-07 ***
lotsize 4.691e-03 3.220e-04 14.569 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3336 on 3021 degrees of freedom
Multiple R-squared: 0.6904, Adjusted R-squared: 0.6901
F-statistic: 2246 on 3 and 3021 DF, p-value: < 2.2e-16
VISÃO GERAL:
Seção de resíduos: avaliação rápida de resíduos. Idealmente, o primeiro e o terceira quartis e Min/Max serão aproximadamente equivalentes em valor absoluto.
Coeficientes: lista os coeficientes estimados juntamente com
testes da hipótese nula de que cada coeficiente é estatisticamente igual
a zero. Est/SE = valor t.
Erro padrão residual: estimativa do desvio padrão constante da distribuição dos resíduos, que por hipótese segue uma distribuição normal.
graus de liberdade: tamanho da amostra - número de coeficientes (3025 - 4)
R-quadrado: porcentagem da variância total de y explicada pelo modelo.
Estatística F: teste geral de que todos os coeficientes (exceto intercepto) são 0.
Todos os valores-p referem-se aos testes de hipótese de que os coeficientes são iguais a zero. Muitos estatísticos e pesquisadores preferem observar os intervalos de confiança.
round(confint(m2) * 100, 4)
2.5 % 97.5 %
(Intercept) 1157.2037 1166.5736
finsqft 0.0461 0.0492
bedroom 2.6363 5.9473
lotsize 0.4060 0.5323
De acordo com o nosso modelo, cada quarto adicional acrescenta, em média, entre 2% e 6% ao valor de uma casa, assumindo que todas as demais variáveis preditoras se mantém constante.
Escreva código para modelar log(totalvalue) como
função de fullbath e finsqft. Chame seu modelo
de m3
Escreve o código para produzir os gráficos de diagnóstico dos resíduos
Como interpretamos o coeficiente estimado da variáveil
fullbath?
Escreva código para simular os dados do modelo, em seguida,
compare com o totalvalue observado. Este parece ser um bom
modelo?
Vamos adicionar hsdistrict ao modelo que acabamos de
ajustar. Estar em um determinado distrito escolar afeta o valor total de
uma casa?
A função table produz uma tabela de frequência
(contagem) de variáveis categóricas (e também de números inteiros!):
table(homes$hsdistrict)
Albemarle Monticello Western Albemarle
1196 983 846
Esses níveis não são números, então como R lida com isso em um modelo linear? A linguagem cria um contraste, que é uma matriz de zeros e uns. Se você tiver um fator com K níveis, terá K-1 colunas.
Neste caso teremos duas colunas: uma para Monticello HS e outra para Western Albemarle HS. Por padrão, R pega qualquer nível que venha primeiro em ordem alfabética e o torna o nível baseline ou referência.
Vejamos o contraste padrão, chamado treatment contrast. Para
fazer isso, convertemos a variável/coluna hsdistrict em un
fator e então usamos a função contrasts().
OBS: Não precisamos fazer isso para adicionar hsdistrict ao nosso modelo! Estamos azendo isso apenas para gerar a “definição” de contraste.
contrasts(factor(homes$hsdistrict))
Monticello Western Albemarle
Albemarle 0 0
Monticello 1 0
Western Albemarle 0 1
Um modelo com hsdistrict terá dois coeficientes:
Monticello e Western Albemarle
Vamos ajustar nosso novo modelo.
m4 <- lm(log(totalvalue) ~ fullbath + finsqft + hsdistrict, data = homes)
summary(m4)
Call:
lm(formula = log(totalvalue) ~ fullbath + finsqft + hsdistrict,
data = homes)
Residuals:
Min 1Q Median 3Q Max
-2.81786 -0.15390 0.00654 0.15716 2.75110
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.158e+01 1.731e-02 669.213 < 2e-16 ***
fullbath 1.617e-01 8.752e-03 18.479 < 2e-16 ***
finsqft 3.911e-04 8.673e-06 45.090 < 2e-16 ***
hsdistrictMonticello -6.748e-02 1.391e-02 -4.849 1.30e-06 ***
hsdistrictWestern Albemarle 9.981e-02 1.472e-02 6.782 1.42e-11 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.322 on 3020 degrees of freedom
Multiple R-squared: 0.7117, Adjusted R-squared: 0.7113
F-statistic: 1864 on 4 and 3020 DF, p-value: < 2.2e-16
Os coeficientes para Monticello e Western Albemarle são em relação a Albemarle HS.
round(coef(m4) * 100, 4)
(Intercept) fullbath
1158.1537 16.1734
finsqft hsdistrictMonticello
0.0391 -6.7479
hsdistrictWestern Albemarle
9.9805
Pelas estimativas, o valor de uma casa em Western Albemarle será cerca de 10% superior ao de uma casa equivalente em Albemarle. Da mesma forma, o valor de uma casa no distrito de Monticello será cerca de 7% inferior ao de uma casa equivalente no distrito de Albemarle.
Escreva o código para modelar log(totalvalue) como
função de fullbath, finsqft e
cooling. Chame seu modelo de m5.
Qual é a interpretação do coeficente de
cooling?
Em nosso modelo acima, que incluía hsdistrict, assumimos
que os efeitos eram aditivos. Por exemplo, não importava em que
distrito escolar a casa ficava, o efeito de banheiro ou
finsqft era o mesmo.
Os efeitos serem aditivos, também implica que o efeito de cada “banheiro completo” adicional era o mesmo, independentemente do tamanho da casa, e vice-versa. Isso pode ser muito simplista.
As interações permitem que os efeitos das variáveis dependam de outras variáveis. Novamente o conhecimento do assunto auxilia na proposição de interações. Como veremos, as interações tornam seu modelo mais flexível, mas mais difícil de entender.
R simplifica a inclusão de interações em modelos. Basta indicar uma interação entre duas variáveis colocando dois pontos (:) entre elas. Abaixo incluímos interações bidirecionais. (Você pode ter interações de três vias e superiores, mas elas são muito difíceis de interpretar.)
m6 <- lm(log(totalvalue) ~ fullbath + finsqft + hsdistrict +
fullbath:finsqft + fullbath:hsdistrict +
finsqft:hsdistrict, data = homes)
summary(m6)
Call:
lm(formula = log(totalvalue) ~ fullbath + finsqft + hsdistrict +
fullbath:finsqft + fullbath:hsdistrict + finsqft:hsdistrict,
data = homes)
Residuals:
Min 1Q Median 3Q Max
-2.7519 -0.1560 -0.0088 0.1461 2.8508
Coefficients:
Estimate Std. Error t value
(Intercept) 1.135e+01 3.640e-02 311.961
fullbath 2.258e-01 1.623e-02 13.907
finsqft 5.885e-04 1.896e-05 31.036
hsdistrictMonticello -2.800e-01 3.820e-02 -7.330
hsdistrictWestern Albemarle -1.186e-01 4.074e-02 -2.911
fullbath:finsqft -6.289e-05 4.308e-06 -14.600
fullbath:hsdistrictMonticello 1.167e-01 2.034e-02 5.739
fullbath:hsdistrictWestern Albemarle 1.190e-01 2.092e-02 5.689
finsqft:hsdistrictMonticello -1.614e-05 2.075e-05 -0.777
finsqft:hsdistrictWestern Albemarle -2.112e-05 2.036e-05 -1.037
Pr(>|t|)
(Intercept) < 2e-16 ***
fullbath < 2e-16 ***
finsqft < 2e-16 ***
hsdistrictMonticello 2.94e-13 ***
hsdistrictWestern Albemarle 0.00363 **
fullbath:finsqft < 2e-16 ***
fullbath:hsdistrictMonticello 1.05e-08 ***
fullbath:hsdistrictWestern Albemarle 1.40e-08 ***
finsqft:hsdistrictMonticello 0.43693
finsqft:hsdistrictWestern Albemarle 0.29967
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3099 on 3015 degrees of freedom
Multiple R-squared: 0.7334, Adjusted R-squared: 0.7326
F-statistic: 921.8 on 9 and 3015 DF, p-value: < 2.2e-16
A interpretação é muito mais difícil. Não podemos interpretar
diretamente os efeitos principais de fullbath,
finsqft ou hsdistrict. Eles interagem. Qual é
o efeito de finsqft? Ele depende de fullbath e
hsdistrict.
As interações são “significativas” ou necessárias? Podemos usar a
função anova para avaliar esta questão. Essa função executa
uma série de testes F parciais. Cada linha abaixo é um teste de
hipótese.
A hipótese nula é que o modelo com este preditor é igual ao modelo sem o preditor. Os testes anova abaixo usam o que é chamado de somas de quadrados do Tipo I. Isso respeita a ordem das variáveis no modelo. Assim:
o primeiro teste compara um modelo com apenas um intercepto a um
modelo com intercepto e com a variável fullbath.
o segundo teste compara um modelo com intercepto e
fullbath com um modelo com intercepto,
fullbath e finsqft.
E assim por diante.
Se a hipótese nula for verdadeira, o valor F deve estar próximo de 1.
anova(m6)
Analysis of Variance Table
Response: log(totalvalue)
Df Sum Sq Mean Sq F value Pr(>F)
fullbath 1 532.64 532.64 5546.784 < 2.2e-16 ***
finsqft 1 227.81 227.81 2372.402 < 2.2e-16 ***
hsdistrict 2 12.54 6.27 65.275 < 2.2e-16 ***
fullbath:finsqft 1 17.76 17.76 184.952 < 2.2e-16 ***
fullbath:hsdistrict 2 5.75 2.88 29.950 1.32e-13 ***
finsqft:hsdistrict 2 0.11 0.06 0.586 0.5566
Residuals 3015 289.52 0.10
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
A interação finsqft:hsdistrict não parece contribuir
muito para o modelo.
O fato de uma interação ser significativa não implica necessariamente que seja relvante ou que valha a pena incluí-la no modelo. Não podemos inferir nada sobre a natureza da interação a partir da tabela da Análise da Variância (ANOVA).
Gráficos de efeitos podem nos ajudar a visualizar e dar sentido aos modelos com interações. Vamos fazer um usando o pacote ggeffects e anlisar o que ele está exibindo.
library(ggeffects)
plot(ggpredict(m6, terms = c("fullbath", "hsdistrict")))
Model has log-transformed response. Back-transforming predictions to
original response scale. Standard errors are still on the log-scale.
# coloca fullbath no eixo x, agrupa por hsdistrict
Qual é o efeito de fullbath? Depende. É mais forte em
Western Albemarle e Monticello. É claro que grande parte da diferença
ocorre em valores extremos de fullbath. As “faixas” ao
redor das linhas representam intervalos com 95% confiança .
O que exatamente foi plotado? Podemos ver usando a funcão
ggpredict sem plot:
ggpredict(m6, terms = c("fullbath", "hsdistrict"))
Model has log-transformed response. Back-transforming predictions to
original response scale. Standard errors are still on the log-scale.
# Predicted values of totalvalue
# hsdistrict = Western Albemarle
fullbath | Predicted | 95% CI
-------------------------------------------
0 | 2.44e+05 | [2.27e+05, 2.62e+05]
1 | 3.02e+05 | [2.89e+05, 3.16e+05]
3 | 4.65e+05 | [4.50e+05, 4.80e+05]
4 | 5.77e+05 | [5.44e+05, 6.11e+05]
5 | 7.15e+05 | [6.56e+05, 7.80e+05]
8 | 1.36e+06 | [1.14e+06, 1.63e+06]
# hsdistrict = Monticello
fullbath | Predicted | 95% CI
-------------------------------------------
0 | 2.09e+05 | [1.96e+05, 2.24e+05]
1 | 2.59e+05 | [2.49e+05, 2.70e+05]
3 | 3.97e+05 | [3.85e+05, 4.09e+05]
4 | 4.91e+05 | [4.65e+05, 5.19e+05]
5 | 6.08e+05 | [5.60e+05, 6.60e+05]
8 | 1.15e+06 | [9.74e+05, 1.36e+06]
# hsdistrict = Albemarle
fullbath | Predicted | 95% CI
-------------------------------------------
0 | 2.86e+05 | [2.68e+05, 3.06e+05]
1 | 3.15e+05 | [3.03e+05, 3.29e+05]
3 | 3.82e+05 | [3.73e+05, 3.92e+05]
4 | 4.21e+05 | [4.01e+05, 4.42e+05]
5 | 4.64e+05 | [4.31e+05, 5.00e+05]
8 | 6.19e+05 | [5.30e+05, 7.23e+05]
Adjusted for:
* finsqft = 2057.60
Not all rows are shown in the output. Use `print(..., n = Inf)` to
show all rows.
ggpredict usou nosso modelo para fazer previsões de
totalvalue para vários valores de fullbath nos
três distritos escolares, mantendo finsqft igual a 1828 (a
mediana de finsqft).
Podemos especificar os valores se quisermos. Por exemplo, crie um
gráfico de efeito para 1 a 5 banheiros e mantenha finsqft
em 2.000:
plot(ggpredict(m6, terms = c("fullbath[1:5]", "hsdistrict"),
condition = c(finsqft = 2000)))
Model has log-transformed response. Back-transforming predictions to
original response scale. Standard errors are still on the log-scale.
E quanto aos efeitos de finsqft e fullbath?
Esta é uma interação de duas variáveis numéricas. A segunda
variável deve servir como variável de agrupamento ao criar um gráfico de
efeito. Abaixo definimos fullbath assumir valores de 2 a 5
e finsqft para assumir valores de 1000 a 4000 em passos de
500.
plot(ggpredict(m6, terms = c("finsqft[1000:4000 by=500]", "fullbath[2:5]")))
Model has log-transformed response. Back-transforming predictions to
original response scale. Standard errors are still on the log-scale.
O efeito de finsqft parece diminuir em função do número
de banheiros completos que uma casa tiver. Mas existem poucas casas
grandes com 2 banheiros completos e, da mesma forma, poucas casas
pequenas com 5 banheiros completos. Embora a interação seja
“significativa” no modelo, é claramente uma interação muito pequena.
Escreva código para modelar log(totalvalue) como
função de fullbath, finsqft,
cooling e a interação entre finsqft e
cooling. Chame seu modelo de m7. A interação é
importante?
Visualize a interação usando a função ggpredict. Use
[1000:4000 by=500] para definir o intervalo de
finsqft no eixo x. Quão notável é essa interação?
Até agora assumimos que a relação entre um preditor e a resposta é linear (por exemplo, para uma mudança de 1 unidade em um preditor, a resposta muda em um valor fixo).
Essa suposição às vezes pode ser simplista e pouco realista. Felizmente, existem maneiras de ajustar efeitos não lineares em um modelo linear.
Aqui está um exemplo rápido de dados não lineares simulados: um polinômio de 2º grau.
x <- seq(from = -10, to = 10, length.out = 100)
set.seed(3)
y <- 1.2 + 2*x + 0.9*x^2 + rnorm(100, mean = 0, sd = 10)
nl_dat <- data.frame(y, x)
plot(y ~ x, nl_dat)
É evidente que um modelo linear não funcionará bem para estes dados. A relação entre x ey não é adequadamente capturada por um modelo linear.
Se quiséssemos tentar “recuperar” os coeficeintes que usamos na
simulação desses dados, poderíamos ajustar um modelo polinomial usando a
função poly() na sintaxe da fórmula:
## Codigo apenas efeitos ilustrativos
nlm1 <- lm(y ~ poly(x, degree = 2, raw = TRUE), data = nl_dat)
No entanto, a abordagem recomendada para ajustar efeitos não lineares é usar splines naturais em vez de polinômios. Splines naturais essencialmente nos permitem ajustar uma série de polinômios cúbicos conectados em nós localizados no intervalo de nossos dados.
A opção mais fácil é usar a função ns() do pacote
splines, que vem instalado com R. ns significa “natural
splines”. O segundo argumento são os graus de liberdade
(df). Pode ser útil pensar em df como o número
de vezes que a linha suave muda de direção.
Frank Harrell afirma em seu livro Regression Model
Strategies que 3 a 5 df é quase sempre suficiente. Seu
conselho básico é alocar mais df para variáveis que você
considera mais importantes.
Vamos ver como funciona com nossos dados simulados.
library(splines)
nlm2 <- lm(y ~ ns(x, df = 2), data = nl_dat)
summary(nlm2)
Call:
lm(formula = y ~ ns(x, df = 2), data = nl_dat)
Residuals:
Min 1Q Median 3Q Max
-21.0764 -6.7423 0.1437 7.4876 18.2382
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.746 2.420 25.92 <2e-16 ***
ns(x, df = 2)1 -77.728 5.433 -14.31 <2e-16 ***
ns(x, df = 2)2 91.189 2.959 30.82 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 8.642 on 97 degrees of freedom
Multiple R-squared: 0.9231, Adjusted R-squared: 0.9215
F-statistic: 582.2 on 2 and 97 DF, p-value: < 2.2e-16
É impossível interpretar o resultado da função
summary(). Assim, vamos visualizar o ajuste com um gráfico
de efeito.
library(ggeffects)
plot(ggpredict(nlm2, terms = "x"), add.data = TRUE)
Data points may overlap. Use the `jitter` argument to add some
amount of random variation to the location of data points and avoid
overplotting.
Vamos voltar aos dados das casas e ajustar um efeito não linear para
finsqft usando um spline natural com df igual
a 5. Abaixo, incluímos também hsdistrict e
lotsize e permitimos que finsqft e
hsdistrict interajam.
nlm3 <- lm(
log(totalvalue) ~ ns(finsqft, 5) + hsdistrict + lotsize +
ns(finsqft, 5):hsdistrict,
data = homes
)
A função anova permite avaliar o efeito não linear e a
interação. Algum tipo de interação entre finsqft e
hsdistrict parece estar presente.
anova(nlm3)
Analysis of Variance Table
Response: log(totalvalue)
Df Sum Sq Mean Sq F value Pr(>F)
ns(finsqft, 5) 5 749.09 149.819 1548.1056 < 2.2e-16 ***
hsdistrict 2 12.56 6.280 64.8896 < 2.2e-16 ***
lotsize 1 28.35 28.354 292.9897 < 2.2e-16 ***
ns(finsqft, 5):hsdistrict 10 5.21 0.521 5.3875 6.092e-08 ***
Residuals 3006 290.91 0.097
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Novamente, é impossível analisar os resultados retornados pela função
summary.
summary(nlm3)
Call:
lm(formula = log(totalvalue) ~ ns(finsqft, 5) + hsdistrict +
lotsize + ns(finsqft, 5):hsdistrict, data = homes)
Residuals:
Min 1Q Median 3Q Max
-3.2767 -0.1459 0.0064 0.1591 2.6353
Coefficients:
Estimate Std. Error t value
(Intercept) 11.5109098 0.2607313 44.149
ns(finsqft, 5)1 1.0921758 0.2464975 4.431
ns(finsqft, 5)2 1.3571870 0.2692905 5.040
ns(finsqft, 5)3 1.7872400 0.1520974 11.751
ns(finsqft, 5)4 3.4685083 0.5580154 6.216
ns(finsqft, 5)5 3.1883957 0.3163353 10.079
hsdistrictMonticello -0.8587210 0.3167315 -2.711
hsdistrictWestern Albemarle -0.1411991 0.3047520 -0.463
lotsize 0.0051240 0.0003031 16.908
ns(finsqft, 5)1:hsdistrictMonticello 0.8620728 0.2980655 2.892
ns(finsqft, 5)2:hsdistrictMonticello 0.7218438 0.3288128 2.195
ns(finsqft, 5)3:hsdistrictMonticello 0.8443943 0.1891643 4.464
ns(finsqft, 5)4:hsdistrictMonticello 1.2923992 0.6818569 1.895
ns(finsqft, 5)5:hsdistrictMonticello -0.0699139 0.3723492 -0.188
ns(finsqft, 5)1:hsdistrictWestern Albemarle 0.2011020 0.2879070 0.698
ns(finsqft, 5)2:hsdistrictWestern Albemarle 0.2743639 0.3155112 0.870
ns(finsqft, 5)3:hsdistrictWestern Albemarle 0.3513783 0.1832677 1.917
ns(finsqft, 5)4:hsdistrictWestern Albemarle 0.3743132 0.6564909 0.570
ns(finsqft, 5)5:hsdistrictWestern Albemarle 0.0911532 0.3526484 0.258
Pr(>|t|)
(Intercept) < 2e-16 ***
ns(finsqft, 5)1 9.72e-06 ***
ns(finsqft, 5)2 4.93e-07 ***
ns(finsqft, 5)3 < 2e-16 ***
ns(finsqft, 5)4 5.81e-10 ***
ns(finsqft, 5)5 < 2e-16 ***
hsdistrictMonticello 0.00674 **
hsdistrictWestern Albemarle 0.64317
lotsize < 2e-16 ***
ns(finsqft, 5)1:hsdistrictMonticello 0.00385 **
ns(finsqft, 5)2:hsdistrictMonticello 0.02822 *
ns(finsqft, 5)3:hsdistrictMonticello 8.35e-06 ***
ns(finsqft, 5)4:hsdistrictMonticello 0.05813 .
ns(finsqft, 5)5:hsdistrictMonticello 0.85107
ns(finsqft, 5)1:hsdistrictWestern Albemarle 0.48492
ns(finsqft, 5)2:hsdistrictWestern Albemarle 0.38460
ns(finsqft, 5)3:hsdistrictWestern Albemarle 0.05530 .
ns(finsqft, 5)4:hsdistrictWestern Albemarle 0.56860
ns(finsqft, 5)5:hsdistrictWestern Albemarle 0.79605
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3111 on 3006 degrees of freedom
Multiple R-squared: 0.7322, Adjusted R-squared: 0.7306
F-statistic: 456.5 on 18 and 3006 DF, p-value: < 2.2e-16
Os gráficos de efeitos são nossa única esperança para compreender este modelo. Abaixo, plotamos o valor total previsto em finqft variando de 1.000 a 3.000, agrupado por hsdistrict.
plot(ggpredict(nlm3, terms = c("finsqft[1000:3000 by=250]", "hsdistrict")))
Model has log-transformed response. Back-transforming predictions to
original response scale. Standard errors are still on the log-scale.
O efeito de finsqft no valor total parece
mais forte em Western Albemarle quando você ultrapassa 1.500 pés
quadrados.
O modelo simula dados semelhantes em distribuição aos dados observados?
sim4 <- simulate(nlm3, nsim = 50)
plot(density(log(homes$totalvalue)))
for(i in 1:50)lines(density(sim4[[i]]), col = "grey80")
Ainda devemos verificar as hipóteses sobre os resíduos do modelo.
plot(nlm3)
As casas 12, 40, 963 e 1810 parecem se destacar. Vamos dar uma olhada.
h <- c(12, 40, 963, 1810)
homes[h,c("totalvalue", "finsqft", "lotsize")]
As casas 12 e 40 têm valor total muito baixo e o modelo superestima seus valores. A casa 963 tem um valor total enorme com 0 acres de tamanho de lote. A casa 1810 ocupa 611 acres e esse valor tem uma grande influência em seu valor ajustado.
cbind(observed = homes$totalvalue[h], fitted = exp(fitted(nlm3)[h]))
observed fitted
12 9600 132267.3
40 16400 434393.2
963 3117300 223502.3
1810 2488800 5491200.9
Escreva código para modelar log(totalvalue) como
função de finsqft com um spline natural com
df = 5, cooling, e a interação de
cooling e finsqft (spline natural com
df = 5). Chame seu modelo de nlm4.
Use a função anova para verificar se a interação
parece necessária. O que você acha?
Crie um gráfico de efeito de finsqft por
resfriamento. Tente [1000:5000 by=250] para o
intervalo de valores de finsqft.
Esta atividade teve como objetivo mostrar a vocês os fundamentos da modelagem linear em R. Esperamos que você tenha uma compreensão melhor de como funciona a modelagem linear.
O que fizemos hoje funciona para resultados numéricos
independentes. Tínhamos uma observação por casa e a nossa resposta
foi “valor total” (totalvalue), um número. Nossos modelos
retornaram o valor total médio esperado, dados vários
preditores. Este é um design bastante simples.
As coisas ficam mais complicadas quando você tem, digamos, respostas binárias ou múltiplas medidas sobre a mesma observação. Uma lista não exaustiva de outros tipos de modelos inclui:
modelos lineares generalizados (para respostas binárias e de contagem)
modelos logit multinomiais (para respostas categóricas)
modelos logit ordenados (para respostas categóricas ordenadas)
modelos lineares longitudinais ou de efeito misto (para respostas com múltiplas medidas sobre o mesmo assunto ou grupos de medidas relacionadas)
modelos de sobrevivência (para respostas que medem o tempo até um evento)
modelos de séries temporais (para respostas que apresentam, digamos, variação sazonal ao longo do tempo)
Referências
Na modelagem anterior, aplicamos a transformação logaritmica em
totalvalue para reduzir a assimetria da distribuição e,
portanto, para ajudar a atender às hipóteses do modelo de regressão
linear clássica.
Lembre-se de que, sem a transformação logarítmica, nossos resíduos eram grandes e assimétricos, o que é uma maneira elegante de dizer que nosso modelo não possuía bom ajuste aos dados. Um bom modelo deve ter resíduos relativamente pequenos com dispersão simétrica.
Uma transformação logarítmica fazia sentido por dois motivos:
A variável resposta totalvalue era estritamente
positiva, tinha um limite superior grande e cobria várias ordens de
grandeza.
as mudanças em totalvalue de acordo com os
preditores foram relativas (multiplicativas) e não absolutas (aditivas),
o que corresponde à escala logarítmica natural.
É importante observar que nem todas as variáveis com distribuição assimétrica precisam ser transformadas quando se trata de modelagem linear. As hipóteses distributivas são feitas em relação aos resíduos não em relação às variáveis resposta e preditores.
No entanto, pode haver momentos em que você precise investigar outras transformações além da logarítmica. Essas transformações alternativas, auase sempre assumem a forma de uma transformação de potência (ou seja, eleva-se a variável a uma potência usando um expoente).
As potências são geralmente simbolizadas pela letra grega lambda (\(\lambda\)). Como uma potência igual 0, temos a transformação logarítmica.
Digamos que a variável seja y. Uma paleta básica de
possíveis transformações de poder inclui:
A função symbox() do pacote car cria uma
avaliação visual de qual potência torna a distribuição razoavelmente
simétrica.
Abaixo, quando a usamos em totalvalue, vemos que a
transformação logarítmica (λ = 0) faz o melhor trabalho em tornar a
distribuição mais simétrica.
car::symbox(homes$totalvalue)
Também podemos usar symbox() em um objeto modelo. Por
exemplo, isso produz essencialmente o mesmo gráfico usando os resíduos
do modelo em vez do valor total. Simplesmente canalize o
modelo para symbox().
lm(totalvalue ~ finsqft + bedroom + lotsize, data = homes) |>
car::symbox()
Uma “busca” pela “melhor” transformação de potência pode ser
realizada com a função powerTransform(), também no pacote
car. A prática usual é converter o resultado para a potência simples
mais próxima listada acima. Por exemplo, podemos canalizar o resultado
do modelo para powerTransform() e ver que a “melhor”
transformação é cerca de 0,16.
lm(totalvalue ~ finsqft + bedroom + lotsize, data = homes) |>
car::powerTransform()
Estimated transformation parameter
Y1
0.1616362
0,16 está próximo de 0, então faz sentido prosseguir com uma transformação logarítmica. Isso simplifica muito a interpretação.